/*=============================================================================
*
*   Copyright (C) 2020 All rights reserved.
*
*   Filename:		test_dirichlet1.cpp
*
*   Author: Wang Zhecheng - wangzhecheng@yeah.net
*
*   Date: 2020-02-20 21:23
*
*   Last Editors: Wang Zhecheng - wangzhecheng@yeah.net
*
*   Last modified:	2020-07-20 16:20
*
*   Description:
*
=============================================================================*/
#define BOOST_TEST_MODULE test_dirichlet

#include "../../include/fem.hpp"
#include "../../include/implicit_solver.hpp"
#include "../../include/incremental_load.hpp"
#include "../../include/model/model.hpp"
#include "../../include/solid.hpp"
#include "../../include/static_equation.hpp"
#include <boost/test/included/unit_test.hpp>

LOGGING_REGISTER()

BOOST_AUTO_TEST_CASE(test_dirichlet) {
  /*  detail, trace, process, info, warning, debug, error }; */
  geoxfem::log::log_init(geoxfem::log::severity_level::error);

  using namespace geoxfem;
  using namespace geoxfem::model;
  using namespace geoxfem::geometry;

  using MethodPolicy   = fem;
  using CouplePolicy   = solid;
  using DimensionType  = plane_strain;
  using EquationPolicy = static_equation;
  using SolvePolicy    = implicit_solver;
  using LoadPolicy     = incremental_load;

  geoxfem::model::Model<MethodPolicy, CouplePolicy, DimensionType,
                        EquationPolicy, SolvePolicy, LoadPolicy, mpl1>
    model;

  BOOST_TEST_MESSAGE("Dirichlet Boundary Check #1.\n"
                     "Square object.\n"
                     "    Left : Dirichlet boundary moving rightwards.\n"
                     "    Others: free boundary");
  std::string input_file_name =
    (boost::unit_test::framework::master_test_suite().argc < 2)
      ? "test/test_bc/dirichlet1.json"
      : boost::unit_test::framework::master_test_suite().argv[1];

  model.read(input_file_name);
  model.buildProblem();
  model.solveProblem();

  MatrixXR correct_ans(12, 2), ans(12, 2);
  correct_ans << 0.0000E+00, -0.8601E-05,  //
    0.0000E+00, -0.4056E-05,               //
    0.0000E+00, 0.0000E+00,                //
    -0.5472E-07, -0.3771E-05,              //
    0.1225E-05, -0.1910E-05,               //
    0.0000E+00, 0.0000E+00,                //
    0.1792E-06, 0.5860E-06,                //
    0.6280E-06, 0.1707E-06,                //
    0.0000E+00, 0.0000E+00,                //
    0.0000E+00, 0.1129E-06,                //
    0.0000E+00, 0.1066E-06,                //
    0.0000E+00, 0.0000E+00;
  for(auto& v : shapesIn(*model.getMesh(), 0))
    ans.row(v->getIndex()) << v->toVertex().getDof("x0")->getCurrValue(),  //
      v->toVertex().getDof("x1")->getCurrValue();
  std::cout << "ans norm:" << correct_ans.norm() << ", " << ans.norm()
            << std::endl;
  BOOST_CHECK_CLOSE(correct_ans.norm(), ans.norm(), 1e-2);

  auto shp1 = model.getMesh()->getShape(2, 5);

	arithmetic::Stress<DimensionType> s_shp1gp1, s_shp1gp2, s_shp1gp3, s_shp1gp4,
    correct_s_shp1gp1, correct_s_shp1gp2, correct_s_shp1gp3, correct_s_shp1gp4;
  s_shp1gp1 = shp1->getGaussPoint(0).getCurrSigma();
  correct_s_shp1gp1 << -0.4e-2, 0.9688e-2, 0,  //
    0.9688e-2, 0.2469e-1, 0,                   //
    0, 0, 0;
  BOOST_CHECK_CLOSE(correct_s_shp1gp1.norm(), s_shp1gp1.norm(), 2e-2);

  s_shp1gp2 = shp1->getGaussPoint(1).getCurrSigma();
  correct_s_shp1gp2 << -0.5281e-1, 0.8267e-2, 0,  //
    0.8267e-2, 0.3775e-02, 0,                     //
    0, 0, 0;
  BOOST_CHECK_CLOSE(correct_s_shp1gp2.norm(), s_shp1gp2.norm(), 1e-2);

  s_shp1gp3 = shp1->getGaussPoint(2).getCurrSigma();
  correct_s_shp1gp3 << -0.4855e-1, 0.3616e-1, 0,  //
    0.3616e-1, 0.1373e-1, 0,                      //
    0, 0, 0;
  BOOST_CHECK_CLOSE(correct_s_shp1gp3.norm(), s_shp1gp3.norm(), 1e-2);

  s_shp1gp4 = shp1->getGaussPoint(3).getCurrSigma();
  correct_s_shp1gp4 << 0.264e-3, 0.3758e-1, 0,  //
    0.3758e-1, 0.3465e-1, 0,                    //
    0, 0, 0;
  BOOST_CHECK_CLOSE(correct_s_shp1gp4.norm(), s_shp1gp4.norm(), 1e-2);
}

